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Abstract 

In this contribution we present a new computational method for coupled bulk-surface 
problems on time-dependent domains. The method is based on a space-time formulation 
using discontinuous piecewise linear elements in time and continuous piecewise linear ele¬ 
ments in space on a fixed background mesh. The domain is represented using a piecewise 
linear level set function on the background mesh and a cut finite element method is used 
to discretize the bulk and surface problems. In the cut finite element method the bilinear 
forms associated with the weak formulation of the problem are directly evaluated on the 
bulk domain and the surface defined by the level set, essentially using the restrictions of the 
piecewise linear functions to the computational domain. In addition a stabilization term is 
added to stabilize convection as well as the resulting algebraic system that is solved in each 
time step. We show in numerical examples that the resulting method is accurate and stable 
and results in well conditioned algebraic systems independent of the position of the interface 
relative to the background mesh. 


1 Introduction 

Problems involving phenomena that take place both on surfaces (or interfaces) and in bulk do¬ 
mains occur in a variety of applications in fluid dynamics and biology. In this paper, we consider 
a coupled bulk-surface problem modeling the evolution of soluble surfactants. A soluble surfac¬ 
tant is dissolved in the bulk fluid but also exists in adsorbed form on the interface separating two 
immiscible fluids. Surfactants have a large influence on the dynamics in multiphase flow systems 
in that they may cause drop-breakup or coalescence due to their ability to reduce the surface 
tension. They were for example used to lower the surface tension of oil droplets in the 2010 
Deepwater Horizon oil spill so that the oil became more soluble in the water. Other examples 
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of applications where the effects of surfactants are important include drug delivery, treatment 
of lung diseases, and polymer blending [1]. 

We consider a coupled system of time-dependent convection-diffusion equations describing 
the concentration of surfactants in the bulk ffuid and on the interface. The interface is moving 
with a given velocity. From a computational point of view, the main challenge is that the 
differential equations are defined on domains that are evolving with time and that these domains 
may undergo strong deformations. 

A common strategy is to let the mesh conform to the time-dependent domain, see, e.^., Eli. 
This technique can be made accurate but requires remeshing and interpolation as an interface 
evolves with time and leads to significant complications when topological changes such as drop- 
breakup or coalescence occur, especially in three space dimensions. Therefore, computational 
methods that allow the interface to be arbitrarily located with respect to a fixed background 
mesh, so called fixed grid methods, have become highly attractive and significant effort has 
been directed to their development, see, e.^., m El El El. In fixed grid methods a strategy for 
solving the bulk Partial Differential Equation (PDE) defined on a domain with the interface as 
boundary is to extend the PDE to the whole computational domain by for example regularized 
characteristic functions, cf. BM- Strategies for solving quantities on evolving surfaces are in 
general developed on the basis of the interface representation technique. In consequence, existing 
fixed grid methods are usually tightly coupled to the interface representation. Techniques to 
represent the interface can be roughly divided into two classes: explicit representation, e._g., by 
marker particles m, and implicit representation, c.g.^ by the level set of a higher dimensional 
function m- Existing methods using implicit representation techniques generally extend the 
surfactant concentration to a region embedding the interface, and instead of a surface PDE, 
a PDE on a higher dimensional domain must be solved for the interfacial surfactant. Several 
methods have been proposed, based on explicit naEamiEsiiisiE] as well as implicit HZlEHl 
[niEniEi] representations. Most work has been done on insoluble surfactants, i.e., surfactants 
that are only present at the interface without surfactant mass transfer between the interface and 
the bulk. 

In this paper, we present a new computational method for solving coupled bulk-surface 
problems on time-dependent domains. The surface PDE is solved on the interface which can 
be arbitrarily located with respect to the fixed background mesh. The method is accurate and 
stable and results in well conditioned linear systems independently of how the interface cuts 
through the background mesh, and the total mass of surfactants is accurately conserved. 

Our strategy is to embed the time-dependent domain where the PDE has to be solved in 
a fixed background grid, equipped with a standard finite element space, and then take the 
restriction of the finite element functions to the time-dependent domain. This idea was first 
proposed for an elliptic problem with a stationary fictitious boundary in [7] and for the Laplace- 
Beltrami operator on a stationary interface in [22]. It has been extended to other equations with 
error analysis, for example the Stokes equations involving two immiscible incompressible fluids 
with different viscosities and with surface tension [23], to PDE:s on time-dependent surfaces in, 
e.g., [SIESIEHIEZ], and to stationary coupled bulk-surface problems with linear coupling terms 
in [281 EH]- These types of methods are referred to as cut finite element methods (CutFEM), 
since the interface cuts through the background grid in an arbitrary fashion. 
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We suggest a CutFEM based on a space-time approach with continuous linear elements in 
space and discontinuous piecewise linear elements in time. The method presented in [25] is for 
solving surface PDEs but is also based on a space-time approach with discontinuous elements in 
time. However, in our approach we add a consistent stabilization term jsniEi] which ensures that 
1) our method leads to linear systems with bounded condition number, 2) the discretization of 
the surface PDE is stable also for convection dominated problems, and 3) the proposed method 
relies only on spatial discretizations of the geometry at quadrature points in time and results 
in the same computations as in the case of stationary problems. In addition, the total mass 
of surfactants is accurately conserved using a Lagrange multiplier. Numerical results indicate 
that the method is optimal order accurate (second order); we have also proven optimal order 
of accuracy for a related stationary coupled bulk-surface problem with a linear coupling term 
in [28]. 

In this paper, we have used the standard level set method to represent the interface, but 
other interface representation techniques can be used as well. We have chosen to concentrate 
on the challenging task of solving the coupled bulk-surface problem on time-dependent domains 
and will throughout the paper assume that the velocity field is given. 

The remainder of the paper is outlined as follows. In Section [^ we formulate the coupled 
bulk surface problem. In Section [^ we introduce a discrete approximation of the interface and 
state our assumptions on the geometry. The computational method for the coupled bulk-surface 
problem modeling soluble surfactants is presented in Section]^ In Sectionwe show numerical 
examples and in Section we summarize our results. 

2 The coupled bulk-surface problem 

2.1 The domain 

Let be an open bounded domain in R^, d = 2,3, with convex polygonal boundary and 
let I = [0,T] be a time interval. We consider two immiscible incompressible fluids that occupy 
time dependent subdomains ^i{t) C il, i = 1,2, with t ^ such that Q = U ^ 2 {t) and 

rii(t) n = 0, and are separated by a smooth interface defined by r{t) = dQi{t) H dQ 2 {t). 

See Fig. [^for an illustration in two dimensions. We assume that F(t) is a simple closed curve 
(2D) or surface (3D) moving with a given divergence free velocity field f3 : I x Q ^ in 
such a way that it does not intersect the boundary dD (F H dD = 0) or itself for any t G /. 
In applications, the velocity field (3 is typically obtained from the incompressible Stokes or the 
incompressible Navier-Stokes equations. For simplicity, we further assume that the surfactant 
is soluble only in the outer fluid phase Di(t). 

2.2 The bulk-surface problem 

The model for soluble surfactants is given by a time-dependent convection-diffusion equation 
on the interface coupled with a time-dependent convection-diffusion equation in the bulk. The 
concentration of surfactants in the bulk and the concentration of surfactants on the surface are 
coupled through a nonlinear term which enters as a source term in the PDE for the interfacial 
surfactant and in a Neumann boundary condition for the bulk concentration. More precisely. 
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Figure 1: Illustration of the domain Q G occupied by two immiscible fluids separated by an 
interface F. Immiscible incompressible fluids with density pi and viscosity pi occupy subdomains 
^i{t) C i = 1,2. 


we consider the following time dependent coupled bulk-surface problem: find : / x ^ R 


and : / X F ^ R such that 

dfUB + /3 • Vxx^ — V • {kB^UB) =0 in / x (2-1) 

71 • kjB^U b — /coupling (2*2) 

—non • ^b^ub = 0 on 512 (2.3) 

dtus + fi • Vus + (divr/3)xx5 - divp {ksVrus) = /coupling on / x T{t) (2.4) 


with initial condition xx^(0,x) = xx^ and xx5(0,x) = xx^ on 12i(0) and F(0). Here V is 

the usual R^ gradient, Vp is the tangent gradient associated with F defined by 


Vp = Pr ■ V, Pr — I — 71 ® n (2-5) 

71 is the unit normal vector to F, outward-directed with respect to 12i, tiqq^ is the outward directed 
unit normal vector on 512, I is the identity matrix, (g) denotes outer product (a (g) b)ij — aibj 
for any two vectors a and 5), and kB and ks are the bulk diffusion and the interfacial diffusion 
coefficients, respectively. The divergence divpi; on F(t) of a vector valued function v is defined 
by 

divpi; = tr(i; (g) Vp) (2.6) 

where (i; (g) Vp)^j = (Vp)ji;^. 

The exchange of surfactants between the interface and the bulk is modeled with the term 
/coupling- We consider, in particular, the Langmuir model where 

/coupling k^UB (xx^ '^S') k^Us (2-7) 
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Here is the maximum surfactant concentration on F and ka and are adsorption and 
desorption coefficients, respectively. Examples of other models are 

/coupling kaUB k^us (Henry) (2*^) 

/coupling = kaUB ( 1 “ ^ “ kae^^^ug (Frumkin) (2.9) 

\ “5 / 

see for example [32]. For a numerical study of different isotherms see [33|. Using the fact that 
V • /3 = 0 we also have the conservation law 

/ usdv-\- / usds^uo for t > 0 (2.10) 

jQiit) Jr{t) 

for the total amount of surfactants on the surface and in the bulk. See the Appendix for the 
formulation of the transport equations for soluble surfactants in non-dimensional form. 


2.3 Weak form 


First note that we can write /coupling? defined in equation (2.7), in the form 


,oo 


/coupling — bsUB 

bs = kd, and bss = ka- 


- bgus — bssUBUs ( 2 - 11 ) 

We assume that bs, bg, and bss are positive 


where bs = kaU% 

constants. For bBg = 0 we obtain the coupling term in the Henry case, equation ( |2.8[ ), with 
bs = ka and bg = ka- 

Let W = H^{yL\{t)) X H^{T{t)). Multiply equation ( |2.1[ ) with a test function bBVB € 
H^{Q,i{t)) and equation \2A\ with a test function bgvg E H^{r{t)). After integration by parts 
and incorporating the boundary conditions (2.2) - (2.3) we obtain the variational problem: find 
u = {ub, ug) € W such that 


bB{dtUB,VB)ni{t)+bg{dtUg, Vg)Y(t)+(l{t, U, v)-{bBgUBUg, bBVB-bgVg)Y(t) = 0 


Here 


a{t, u, v) = bBttBit, UB, vb) + bgag{t, ug, vg) + aBs{t, u, v) 


Vu = {vb, Vg) € W 
( 2 . 12 ) 

(2.13) 


with 


aB{t,UB,VB) = (/3 • V'Us,us)Qj(t) + (/cfiVus, Vub)q^(j) 

< ag{t, ug, vg) = {(3 ■ Vug, vg)r(t) + {{diYrP)ug, vg)r(t) + {kgVrug, Vrus)^^^) (2.14) 

^aBs{t,U,v) = [bBUB - bgUg,bBVB - bgVg)Yit) 


Note that aBg{t,u,v) is the linear part of the coupling term /coupling- 
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3 The level set representation of the interface 


For t ^ let U{r{t)) be an open neighborhood of r{t) such that for each x G U{r{t)) there is 
a uniquely determined closest point in T{t). We let p(t^x) : ^ R be the signed distance 

function. The exterior unit normal n = n(t, x) is the spatial gradient of the signed distance 
function, n{t^x) — Vp{t^x)^ for x G r(t). 

Given a vector field /3 the evolution of the surface is governed by the following problem 
for the level set function: find p : I x Q such that 

pt + (3 ■Vp = 0, p{0,x) = po (3.1) 


Let ICo^h be a quasiuniform partition of ft into shape regular triangles for d = 2 and tetrahedra 
for d = 3 of diameter h. We approximate the level set function p by G where V^o,/i/2 is 

the space of piecewise linear continuous functions defined on the mesh /Co,/i/2 obtained by refining 
uniformly once. Th{t) is the zero level set of ph{t^x). We consider a continuous piecewise 
linear approximation F/^ of F such that F/^ H iF is a linear segment for d = 2 and is a subset of 
a hyperplane in for each K G /Co,/i/2* assume that for every t G /, F/j,(t) C U{T{t)) and 
that the following approximation assumptions hold: 


iL^irh) 


<h^ 


(3.2) 


and 

- nh\\L^(Th) Sh ( 3 - 3 ) 

for all t ^ I. Here < denotes less or equal up to a positive constant, Uh denotes the piecewise 
constant exterior unit normal toT^ with respect to fti and is the extension of the exact normal 
to F/i by the closest point mapping. These assumptions are consistent with the piecewise linear 
nature of the discrete interface. We define fth,i as the domain enclosed by F/^ U dft. 

Let 0 = to<ti<-**<tw = Tbea partition of the time interval / = [0, T] into time steps 
In — {tn-i^tn] of length kn — tn — tn -1 for n = 1, 2,..., A^. To solve the advection equation (3.1) 
we use the Crank-Nicolson scheme in time and piecewise linear continuous finite elements with 
streamline diffusion stabilization in space. We obtain the method: find p^ G y^^hl2 such that, 
for n = 1 , 2 ,..., A^, 


^71 1 ^71 1 

+ C^ + -r • tsd(3^ ■ = 

tZji L Kji L 

= • wrt + rsDl3^ • v< G Vo^u/2 (3.4) 

K>yi Z 

where the streamline diffusion parameter tsd = To keep the level set 

function a signed distance function, the reinitialization equation, equation (15) in [3l], can be 
solved, e.p., in the same way as we did with the advection equation in ( |3.4| ). 


4 The space-time cut finite element method 

In a space-time formulation the variational formulation is written over the space-time domain 
which is divided into space-time slabs corresponding to each time step. We employ piecewise 
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linear trial and test functions that are continuous in space but allowed to be discontinuous from 
one space-time slab to another. For a given time tn (where n is the time step number) we thus 
have two distinct solutions, at times := lime^o^n ± ^5 see, e.g. [35]. Using a suitable weak 
enforcement of continuity at tn the discrete equations can then be solved one space-time slab 
at a time. The discontinuous Galerkin method in time can be compared with implicit finite 
difference methods and good stability is expected (the method is related to the first subdiagonal 
Fade approximation, cf. Thomee [36|). In space we will use a CutFEM, which means that we 
will use restrictions of standard continuous finite element functions defined on the background 
grid to our time dependent domain. Since the interface can cut through the fixed background 
grid arbitrarily there is a lack of shape regularity which results in ill-conditioned system matrices. 
To avoid this problem, we add stabilization terms similar to [28]. The same stabilization terms 
also stabilize the proposed finite element method in case of convection dominated problems. 


4.1 The mesh and finite element spaces 

We define the following sets of elements 


ICB,h{t) = {Ke /Co,/. : K n f2/.,i(/) ^ 0}, 
and the corresponding sets 

■^B,h = U U 

t^In KeJCB,hit) 


iCs,k{t) = {Ke K^.k ■■ K n ri(t) 0} (4.1) 


^1. = u u 

tGin Ke)Cs,h{t) 


(4.2) 


For an illustration of the sets in two dimensions see Fig. 

Let Pi{In) be the space of polynomials of order less or equal to 1 on In and let Vb,/i be the 
space of piecewise linear continuous functions defined on JCo^h- On each of the space-time slabs 
— In'X- AAg ^ and — In x we define the spaces 


= Pi (4) ® Uo,/.k 


S,h 


and we let 




Functions in take the form 


/ \ / \ I ^ ^Tl —1 ^ —1 i 

v[t,x) = [VB^VS) = [VB,0 + VB,1 - 7 - ^VS,0 + VS,1 - 7 - 

\ J 

where t ^ In and vbj and vsj^ j = 0,1 can be written as 

Nb Ns 

VB,j = = Y \^s,h 

i=l i=l 


(4.3) 

(4.4) 


(4.5) 


(4.6) 


Here G R are coefficients, g:)i{x) is the standard nodal basis function associated with 

mesh vertex i, Nb and Ns are the number of nodes in Mq ^ and in respectively. 
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Figure 2: Illustration of the sets introduced in Section 4.1. In both figures the blue and the 
red curves show the position of the interface at the endpoints t — tn-i and t of the time 
interval In — Left: the shaded domain shows AAg ^ and edges in ^ are marked 

with a thick line. Right: the shaded domain shows and edges in J^s,h are marked with a 
thick line. 


4.2 The finite element method 

Given and Hq (see equation ( 2.10[ )) find Uh — {ub^us) G and A G R, such that 

Vh) + Jh{uh, Vh) + A (^(1, + (1, + M ((«b, l)oh,i(tn) + («5, l)r,.(t„) 


= /iMo, 'ivh € W^, e R 


(4.7) 


Here 


^hi'>^^v)= bBidtUB,VB)nf^^^(t)dt+ bsidtus,vs)rf,(t)dt + ahit,u,v)dt 
J JJ 

- / bBs{uBUs,bBVB -bsvs)ri,(t)dt 
Jin 

+ bB{[uB],VB{tn-i,x))Q^^^(^tn-i) + f^s{[us],vs{t+_^,x))r,^(tn-i) (4-8) 


with 

and 


ah{t, u, v) = bBaB,h{t, ub, vb) + bsas,h{t, us, vs) + aBS,hit, u, v) 


(4.9) 


0-B,h{t, UB, Vb) = (/3 • VUB, VB)Qy^^^(t) + {kB'^UB, 

as,h{^, US, vs) = (/3 • '^us,vs)r^{t) + ((divr;,/3)«5,t'5)r,.(t) + {ks'^rHUs,'^rHVs)rf,^(t) 

aBS,h{t,u,v) = {bsUB — bsus,bBVB — bsvs)rh{t) 


(4.10) 
























































































where Vr^, = Ph ’ V and — I — rih® Next 

Jh{u,v)=[ jh{u,v)dt (4.11) 

Jin 

where jh{u^v) is a stabilizing term of the form 


jh{v,w) = TBhjB{vB,WB) + Tsjs{vs,ws) 


(4.12) 


iBiTS are positive parameters and, letting [a:]|F denote the jump of x over the face F, 


jB{vB,WB) — ^ {[np ■'VvB],[nF ■'S/wb])f (4-13) 

js{vs,ws) = ^ ([np • Vns'], [np • Vwf])^ (4-14) 


with Ps,h fhe set of internal faces (i.e. faces with two neighbors) in and pB,h fhe set of 

see Fig. Finally, note 


faces that are internal in ^ and also belong to an element in 
that we use a Lagrange multiplier to impose the condition (2.10). 


Remark 4.1. The stabilization terms js and js that appear in the method are eonsistent and 
are needed to eontrol the eondition number of the system matrix so that the resulting algebraie 
system is well eonditioned independently of the position of the interfaee relative to the mesh. 
The stabilization term js also ensures a stable diseretization of the surfaee PDE in ease the 
problem is eonveetion dominated m- However, one may need to apply the stabilization term 
jB on all faees that are internal in J^b h guarantee a stable diseretization of the bulk PDE 
in ease the problem is eonveetion dominated. One may also apply the bulk stabilization only on 
faees in Tb^h ond add other stabilization methods like for example SUPG HI when the problem 
is eonveetion dominated. 


Remark 4.2. The proposed CutEEM method with eontinuous linear elements in spaee and 
diseontinuous pieeewise linear elements in time is seeond order aeeurate both in spaee and time. 


Remark 4.3. We may also eonsider higher order diseontinuous Galerkin method in time based 
on polynomials of order p. The optimal order of eonvergenee for a parabolie problem is p + 1 
inside the intervals In and 2p + 1 in the nodes tn, see m for further details. Note, however, 
that in order to aehieve higher order eonvergenee in time both the transport equation for 

the distanee funetion p and the bulk-surfaee problem (2. 1-2.4) for the eoneentrations ub and us 
must be diseretized with a higher order method. 


Remark 4.4. We may eonsider the same method without the Lagrange multipliers A and p. 
This method is also of optimal eonvergenee order but the eonservation of the total mass of 
surfaetants is lost. See also Example 1 in Seetion\^ Strong imposition of the eonservation law 
using Lagrange multipliers essentially eompensates for numerieal errors, sueh as the error in the 
area of the surfaee and in the volume of the bulk domain, during eaeh time step. 
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4.3 Implementation 
4.3.1 Newton’s method 


Since the bulk and surface surfactant forms are coupled through a nonlinear term, see (2.7), 


the proposed method ( |47l ) leads to a nonlinear system of equations in each time step, which we 
solve using Newton’s method. To formulate Newton’s method we define the residual F and the 
Jacobian DF as follows 

F{u,\)= / hB{dtUB,VB)Q^,^{t)dt+ / hs{dtus,vs)Th(t)dt+ / ah{t,u,v)dt 
J J J 

- / bBsiuBUs,bBVB - bsvs)rf,{t)dt 
din 

+ bB{[uB],VB{tt-vx))n^^(tn-i) + bs{[us],vs{t^_i,x))r^(tn-i) + / jh{u,v)dt 

Jin 


+ 


A ((1, VB)af,^^{tn) + ^shhitn)) + + (us, l)rh(tn)) “ (4-15) 


, v) dt 


DF{u,\){w,\) = / bBidtWB,VB)nf^^^(t)dt+ / bsidtws,vs)rf,(t) dt + / ahit,w, 

JJ J 

- / hBsiwBUg, bsVB — bsvs)rH{t) dt- bssiuBWs, bsVB - bsvs)r,^{t) 

JJ 

bB{wB,VB{t^_i,x))a^^(tn-i) + bs{ws,vs{tt-i,x))r^{tn-i) + / jh{w,v) 

din 

A ((1, VBh^^^itn) + (1> vs)rh{tn)) + [{WB, + (ws, l)rh(tn)) 


+ 


dt 

dt 


(4.16) 


With this notation the nonlinear problem resulting from (4.7) takes the form: find Uh G Wi 


and A G R such that F{uh^ A) = 0, and the corresponding Newton’s method reads: 

1. Choose initial guesses Uh^o and Aq 

2. while 11 (re, A)|| > tol 

• Solve: DF{uh,o,Xo){wA) = F{uh,oAo) 

• Update Uh^o- = Uh^o — w and Aq: Aq = Aq — A 

For t ^ In we choose the initial guess to be the solution at t~_^ i.e. 

4.3.2 Assembly of the bilinear forms using quadrature in time 

Recall from equation (4.5) that a function v G WJ^ can be written as 

v{t,x) = {vB,vs) = [vb^ + vb / + {Fll) 


10 







where t ^ In and vb,o and vb,i are functions in (fhe space of restrictions to of 

functions in Vo^h) 2 tnd vs^o and vs^i are functions in (fhe space of restrictions to J^Bh 

of functions in Vb,/^). Thus, for example the first term in DF{u^ A)(r(;, A) can be written as 


f bB{dtWB,VB)nf, :^{t)dt = [ dt + [ 

I In I In ^ I In 


t tn—1 / 

bs —p— 




dt 


(4.18) 

Our approach is to first use a quadrature rule in time and for each of the quadrature points 
compute the integrals in space. Since the geometry changes in time the contributions to the 
Jacobian and the residual are computed in each quadrature point for the current geometry 
at that point in time. Using a quadrature formula in time with quadrature weights and 
quadrature points g = 1,..., where riq is the number of quadrature points, the first term 
in DF{u^ A)(r(;, A) is approximated by 


/ bsidtiUB, vb)q^ i(t) dt 
Jin 


E 

q=l 




bs, 


u;-„' — {wB,i- 




F — tn-l 

-p- 

q=l ^ 


+ E 


iwB,l,VB,l)af^,(tq) 

(4.19) 


The other terms in F{u^ A) and DF{u^ A) (it;. A) are treated in the same way. In the numerical ex¬ 
amples in Section we use both the trapezoidal rule and Simpson’s rule for the time integration. 
Choosing the trapezoidal rule for the time integration corresponds to choosing = 2, quadra¬ 
ture points ti = tn-i and ^2 = tn-> and quadrature weights 00^ — Choosing Simpson’s 

quadrature rule corresponds to choosing riq = 3, quadrature points = tn-i^ ^2 = 

For higher order discontinuous 


fU 


Lq 

tn^ and quadrature weights ^ and 


Galerkin methods in time, see Remark [4.3| it is convenient to use the Gauss-Lobatto quadrature 
rule, which includes the endpoints of the interval as quadrature points and is exact for polyno¬ 
mials of order 2nq — 3 and thus have an error term of the form 0{kn ^ )• The optimal local 

order in a discontinuous Galerkin method in time is 2pF 2, where p is the order of polynomials, 
and therefore it is natural to choose riq — p + 2. When using a quadrature rule that includes the 
endpoints some computations can be reused when passing from one space-time slab to another. 

For each of the quadrature points in the time interval In we compute the discrete surface 
defined as the zero level set of the approximate distance function Ph{tq). The intersection 
Th{tq) n K is then planar, since Ph{tq) is piecewise linear, and we can therefore easily compute 
the contribution to the stiffness matrix using a quadrature rule for a line segment in two dimen¬ 
sions and triangles in three dimensions. The contribution from integration on I}h^i{tq) H iF is 
divided into contributions on one or several triangles in two dimensions and tetrahedra in three 
dimensions depending on how the interface cuts element K. 

Note that when we use Simpson’s quadrature rule we need to compute Ph{tn)^ Ph{tn+i)^ 
and Ph{tn+i/2) with = F-i+tn ^ Therefore we use half the time step size i.e. k/2 when 

evolving the interface (solving the advection equation for the level set function) while the coupled 
bulk-surface problem is solved with time step size k. 

Finally, we obtain a 2{Nb + Ns) + 1 x 2{Nb + Ns) + 1 linear system of equations 


DF{uhp, Xo){w, A) = F{uhp, Ao) 


(4.20) 
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for A G R and 


w = 


WB,0 \ 

wsfl 

Wb,1 

\ ^ 5,1 / 


We use a direct solver to solve this linear system of equations. 


4.3.3 Different models of the bulk-surface coupling 


In the proposed method it is straightforward to account for other forms of the coupling term 
/coupling- We now briefly explain how the proposed method should be modified when other 
examples for /coupling are used. 

The variational formulation and thus F{uh^ A) contains the term jj (/coupling? dt. 

The linear part of /coupling is contained in ass^h and the fourth term in F{uh^ A) (equation ( |4.15[ )) 
is the nonlinear part of /coupling- Thus, these two terms and consequently the Jacobian DF 
changes if the coupling term changes. 


For example, in the Henry case (2.8) ass^h is as before but the fourth term in F{uh^X) 
vanishes since hss — ^ and the problem is linear. In the Frumkin case /coupling is of the form 


/coupling — b^UB bgC b^S'^B'^S 


(4.21) 


Hence ass^h = {bBUB^bBVB — bsvs)rh{t) and the fourth term in equation (4.15) is replaced by 



us + bBSUBUS, bBVB “ bsVs)Th{t) 


and the term 



{bs{Ae^'^^us + e^^^)ws, bsVB - bsvs)rh{t) dt 


has to be added to equation (4.16). 


(4.22) 


(4.23) 


Remark 4.5. Another approach to assembly of the discrete problem, used for example in is 
to explicitly construct three dimensional triangulations, in case of two space dimensions, or four 
dimensional triangulations, in case of three space dimensions, of the space-time subdomains. 
However, in comparison our approach is very simple to implement since it only relies on spatial 
discretizations of the geometry at the quadrature points in time, which is the same computation 
as in the case of a stationary problem. 


Remark 4.6. The implementation of the method is straightforward also in three spatial dimen¬ 
sions, since the use of quadrature in time essentially reduces the geometric computations to the 
corresponding stationary problem. A detailed study of a linear coupled bulk-surface problem, in¬ 
cluding an implementation in 3D, estimates of the error, and estimates of the condition number, 
is presented in f^ . Due to the stabilization terms the condition number is always under control 
enabling the use of efficient linear algebra solvers. 
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Remark 4.7. If a coupled bulk-surface problem is solved where the bulk and surface PDEs, 
equation ( |2.lD and ( |2.4[ )^ have nonzero right handsides fsit^x) and fs{t^x) one has to add the 
terms fj bB{fB^'^B)Q^ ^{t) dt and Jj bs{fs^'^s)rh(t)d^t to F{u^X). This is the case in Example 1 
in Section\^ 

5 Numerical examples 

In all the computations in this section we use a uniform underlying mesh /Co,/i consisting of 
triangles of size h and a constant time step size of the form k = Ch. The stabilization constants 
tb and rs in the stabilization term jh are 10“^. 

We consider examples from jnmsiEi and we also formulate one example for which we know 
the exact solution. In the examples from HSIE] the coupled bulk-surface problem is formulated 
in non-dimensional form. This leads to some minor changes in the weak formulation, see the 
Appendix. 

In this section we show that the proposed method is second order accurate. We study the 
convergence at time t = 0.5. In case the exact solution is known we measure the order of 
convergence by studying ||(nB, exact - 1(0.5) and ||(ns, exact - ^^ 5 , 2 / 1 )Hr,. (0.5) for different 

mesh sizes h. When the exact solution is not known we measure the order of convergence by 
using consecutive refinements of the underlying mesh and study \\{uB^h ~ '^5,2/1)||o/, 1(0.5) 

— '^S',2/i)||r/,(o.5)- This is also how the convergence is studied in [T6l [9]. 

Note that in the computation of \\{us^h — '^5, 2/1) Hr/, (0.5) concentration needs to 

be evaluated at quadrature points on r/i(0.5) which do not have to lie on r2/i(0.5). Thus, we 
need to extend the concentration us^2h{0.5) out from r2/i(0.5) to Th^O.b). There are different 
ways of doing this extension. However, analogously to our analysis in [28], for each quadrature 
point on Th{0-5) (used in the computation of \\{us^h ~ '^S',2/i)||r/,(o.5)) we find the closest point 
on r2/i(0.5) and evaluate us^2h that point. We use the same approach in the computation of 
\\{uB,h - «B,2fe)|b,,,i(0.5)- 

We will also show in this section that the condition number of the algebraic system of 
equations is bounded independently of how the interface cuts the underlying mesh and that the 
total mass of surfactants can be conserved accurately. 

5.1 Example 1 

To study the convergence of our numerical method we now consider a coupled bulk-surface 
problem where the exact solution is given. The interface is initially a circle centered at (0.5, 0.22) 
with radius ro = 0.17 and the velocity field /3 = (7r(0.5 — y)^Ti{x — 0.5)). The interface moving 
with this velocity field is at time t a circle with radius ro = 0.17 centered at 

Xc — 0.5 + 0.28sin(7rt) 

yc — 0.5 — 0.28cos(7rt) (5-1) 
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Figure 3: Results for Example 1. Position of the interface and the bulk concentration on the 
moving interface at time t=0.5, 1.2, 2 for mesh size h — 1/40 = 0.025 and time step size k = 0.5/i. 


We choose the computational domain to be [0, 1] x [0, 1]. The bulk diffusion and the interfacial 
diffusion coefficients are ks = 0.01 and ks = 1. The exact solution is 


ub = 0.5 + 0.4cos(7rx) cos^ny) cos(27rt) 

ub -\- ^ sin(7rx) cos(7ry) cos(27rt)ni + ^ cos(7rx) sin(7r^) cos(27rt)n2 


us = 


1.5 + 0.4cos(7rx) cos^ny) cos(27rt) 


(5.2) 


where 


ni = 


{x - Xc) 


n2 


Viy - vcY + {x- Xc)^ 

_ (y - vc) _ 

\/(y - ycY + {x- xcY 


(5.3) 


and Xc and yc are given in equation The function ub satisfies the interface and boundary 


conditions (2.2) and (2.3) but the bulk and surface PDEs (2.1) and (2.4) are satisfied with right 


hand sides /b and fs^ respectively. Eor the implementation, see Remark 4.7 


The bulk and interfacial concentrations on the moving interface at time t = 0.5,1.2, 2 are 
shown in Fig. [^and[^ respectively. The mesh size /i = 1/40 and the time step size k = 0.5/i. 

In this example we compare the errors using the trapezoidal rule and Simpson’s rule for the 


time integration. In Fig.^andp 
versus mesh size h using me dr 


we show the errors \\{uB-UB,h)\\af,^i(o. 5 ) and ||(^^S'-^^5',fe)||r;,(o.5) 
brent time integration schemes. The time step size is A; = 0.5/i. 
We show results both with and without prescribing the total mass UBdv-\- usds. In this 
example, the total mass is time dependent and not conserved. However, at the end of each time 
interval we can compute (since the exact solution is known) and prescribe the total mass. We 
observe the expected second order convergence of \\{ub — UB,h)\\Qh i(o.5) Wi'^S — ||r/,(o.5) 

in the L? norm both using the trapezoidal rule and Simpson’s rule with and without prescribing 
the total mass. For the bulk concentration the errors with and without prescribing the total 
mass coincide and hence we only show the results without prescribing the total mass in Fig. 
Also, the different time integration schemes give similar results because the error in the space 
discretization dominates. For the interfacial concentration we see in Fig that prescribing the 
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Figure 4: Results for Example 1. Position of the interface and the surface concentration on the 
moving interface at time t=0.5, 1.2, 2 for mesh size h — 1/40 = 0.025. 



Figure 5: Results for Example 1. The error \\{ub — ||o/, i(o.5) measured in the L? norm 

versus mesh size h. The time step size k = 0.5/i. Results using the trapezoidal rule (stars) and 
the Simpson’s rule (circles) for the time integration are shown. The total mass is not prescribed. 
The error when the total mass is prescribed coincide with the shown results. The dashed line is 
proportional to 
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Figure 6: Results for Example 1. The error \\{us —||r/,(o.5 ) measured in the norm versus 
mesh size h. The time step size is = 0.5/i. Squares and stars represent the trapezoidal rule 
with and without prescribing the total mass, respectively. Diamonds and circles represent the 
Simpson’s rule with and without prescribing the total mass, respectively. The dashed line is 
proportional to h^. 


total mass gives smaller errors than not prescribing. We also see that Simpson’s rule gives 
smaller errors than the trapezoidal rule. Therefore, in all the other examples we use Simpson’s 
rule and prescribe the total mass. 

In Fig. [^andj^we show the convergence of the bulk and interfacial concentrations at time 
t — 0.5 keeping the mesh size fixed but varying the time step size. We show results both using 
the trapezoidal rule (stars) and Simpson’s rule (circles). Using the trapezoidal rule, represented 
by stars in the figures, we expect to see second order convergence in time. For the bulk concen¬ 
tration the convergence is around second order before the space discretization dominates. For 
the interfacial concentration the convergence is faster for large time step sizes but as the time 
step decreases the error becomes around second order before the error in the space discretization 
dominates. For the Simpson’s rule we observe the same behavior but now the order of conver¬ 
gence seems to be third order instead of second order. Note that the time t — 0.5 at which the 
errors are measured is a nodepoint in the time discretization and since we use linear polynomials 
in the DG method the optimal order of convergence in time is third order, see Remark |4.3[ 


5.2 Example 2 


Next we consider Example 6 of m, a surface convection-diffusion problem with a moving 
interface modeling insoluble surfactants. Thus, surfactants only exist on the interface. Initially 
the interface F is a circle centered at the origin with radius tq = 1 and the velocity field /3 = 


(^^^^,0). The initial interfacial surfactant concentration us{0^ x^y) = ylrQ + 2. The interfacial 
diffusion coefficient ks = I- The computational domain is chosen as D = [—2, 6.4] X [-2, 2] 
and we use 148 gridpoints in the x-direction and 71 gridpoints in the y-direction which yields 
a mesh size h ^ 0.06. The time step size /c = /i/8 as in [18]. The surfactant concentration on 
the moving interface at times t = 0,1,2 is shown in Fig[^ In Fig. W the relative error in the 
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Figure 7: Results for Example 1. The error \\{ub — llo/, 1 ( 0 . 5 ) measured in the norm 
versus time step size k for the mesh size h — 1/160 = 0.00625. Stars represent the trapezoidal 
rule and circles the Simpson’s rule, respectively. The dashed line is proportional to h^. The 
dashed dotted line is proportional to 



Figure 8: Results for Example 1. The error \\{us ~||r/,(o.5) measured in the norm versus 
time step size k for the mesh size h = 1/160 = 0.00625. Stars represent the trapezoidal rule 
and circles the Simpson’s rule, respectively. The dashed line is proportional to h^. The dashed 
dotted line is proportional to 
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Figure 9: Results for Example 2. Position of the interface and the surfactant concentration on 
the moving interface at time t=0, 1, 2. 


total surfactant mass versus time is shown and we see that the error is of the order of machine 
epsilon. In Fig. we show the relative change of the area enclosed by the interface. We observe 
a change in the area by less than 0.005% at time t = 2 for /i ps 0.06. In the method presented 
in [18] which is based on the standard level set method a surfactant mass loss of 4-5 % and a 
change in the area by 1% for h = 0.04 is observed, see Figs. 8 and 9 in [18]. In [18] the PDE 
governing the evolution of the interfacial surfactant concentration is extended off the interface 
to other level sets in a neighborhood of the interface (the zero level set). 


5.3 Example 3 


We consider the same example as in Section 4.1 of US]. Initially the interface F is a circle 
centered in (0,1) with radius ro = 0.5 and the velocity field /3 = (—l + y^O). The computational 
domain is chosen as = [—1 ,1] X [0,2]. The non-dimensional numbers are Pcs = 10, Pe = 1, 
Bi = 1, cr = 1, and Da = 0.2. The initial surface and bulk surfactant concentrations are 
us{0^x^y) = 0.4 and usiO^x^y) = 2/3, respectively. 

The bulk and interfacial surfactant concentrations on the moving interface at time t = 0.5 


are shown in Fig. 12 for mesh size h = 2/50 and time step size k — 0.625/i. We see in Fig. 13 
that the relative error in the total surfactant mass is of the order of machine epsilon and hence 
the total surfactant mass is accurately conserved. In the method presented in d which is based 
on the segment projection method a surfactant mass loss of around 0.001 % is observed for the 
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Figure 10: Results for Example 2. The relative error in the total surfactant mass, 
versus time. 


Ju(t) ^-47rro 

47rro 



Figure 11: Results for Example 2. The relative change of the area enclosed by the interface 
versus time. 
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Figure 12: Results for Example 3. Position of the interface and surfactant concentrations in the 
bulk and on the moving interface at time t = 0.5 for mesh size h — 2/50 = 0.04 and time step 
size k = 0.625/1. 


same bulk mesh as we have used but 2.25 times finer mesh for the interfacial concentration, see 
Fig. 8 in fT6] . 

we show the convergence of \\{uB,h - ^5,2/i)||o;,,i(o.5) and \\{us^h - 


In Fig. 14 and Fig. 
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us,2h)\\rh{o.5)^ respectively. We observe second order convergence both in the norm and the 
norm. We use the same mesh sizes as used in m for the bulk concentration. The method 
in US] is also second order. It uses Strang splitting to both split the coupled bulk-surface problem 
and to split the advection and the diffusion parts of the convection diffusion equations. In the 
discretization of the bulk PDE regular finite difference stencils are used which require points 
outside of the domain fli. These values are determined such that the boundary conditions are 
enforced using an embedded boundary method. Errors measured in the norm reported in 
Fig. 8 of [16] are smaller than we observe. This could be explained by the fact that in [16] the 
interfacial surfactant is discretized on 2.25 times finer meshes than the meshes we use and the 
approximation of the interface in m is more accurate than the approximation we have used. 
However, in comparison the presented method is very simple to implement both in two and 
three space dimensions. 


5.4 Example 4 

We consider here the same coupled bulk-surface problem as in Section 5.3 of [9]. The initial 
interface is a circle with radius ro = 0.3 centered in (0.1, 0) and the velocity field is given by 




—-(1 -h cos(7rx)) sin(7r^), 2 cos(7ry)) sin(7rx^ 


(5.4) 


The computational domain is chosen as = [— 1,1] X [ — 1,1]. The non-dimensional numbers 
are set to Pe = Pes = 100 and Bi = a — Da = 1. The initial surface and bulk surfactant 
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Figure 13: Results for Example 3. The relative error in the total surfactant mass versus time 
for mesh size h — 2/50 = 0.04 and time step size k = 0.625/i. 



Figure 14: Results for Example 3. The error \\{ub — llo/, i(o.5) measured in the norm 
(circles) and the norm (stars) versus mesh size h. The time step size k = 0.625/i. The dashed 
line is proportional to h^. 
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Figure 15: Results for Example 3. The error \\{us^h ~ '^S',2/i)||r/,(o.5) measured in the norm 
(circles) and the norm (stars) versus mesh size h. The time step size k — 0.625/i. The dashed 
line is proportional to h^. 


concentrations are us{0^x^y) = 0 and 


UB{0,x,y) = 


0.5(1 — x^)^ if r > 1.5ro 


0.5(1 — x^yw{r) 
0 otherwise 


if ro < r < 1.5ro 


with r = \/{x- xq)^ + {y - 2 / 0 ) 


^ and 


w{r) ~ 2 ~ 


(r - ro)7T 
0.5ro 


(5.5) 


The bulk and surface surfactant concentrations on the moving interface at times t = 0.5,1,1.5, 2 


are shown for mesh size /i = 2/64 = 0.03125 in Fig. 16 and Fig. 17, respectively. The time step 
size is /c = h/8. We see in Fig. that the proposed method accurately conserves the total 
surfactant mass. The method in [9] extends the bulk equation from to the whole domain by 
a regularized indicator function, therefore there is a mass leakage to the domain ^2 of tho order 
of the regularization parameter. 

we show the convergence of \\{uB^h ~ '^ 5 , 2 / 1 )llo/, 1 ( 0 . 5 ) (represented by circles) and 
ho.5) (represented by stars) in the norm. We observe second order conver- 


In Fig. 
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\\{us,h - '“5,2/1.)||^^^ 

gence which is optimal since we use linear elements in space. For the two coarsest meshes the 
error \\{uB^h ~'^B,2h)\\nh 1 ( 0 .5) reported in [9] are slightly smaller then the errors we observe how¬ 
ever since the method proposed in [9] is only first order accurate the errors we observe decrease 
faster and for the two finest meshes we obtain smaller errors. However, for the mesh sizes shown 
in the figure the error \\{us^h ~ '^S', 2 / 1 )||r/,(o.5) reported in [9] is smaller than the error we obtain. 
This could be explained by the fact that the interface approximation is more accurate in [9] with 
a set of Lagrangian markers. 

We see in Fig. where the condition number is shown as a function of time that as the 
interface evolves the condition number is bounded, independently of how the interface cuts 
through the mesh. 
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Figure 16: Results for Example 4. Position of the interface and the bulk concentration at time 
t — 0.5,1,1.5, 2 for mesh size /i = 2/64 = 0.03125 and time step size k = h/8. 
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Figure 17: Results for Example 4. Position of the interface and the surface concentration at 
time t — 0.5,1,1.5, 2 for mesh size h — 2/64 = 0.03125 and time step size k — h/^. 



Figure 18: Results for Example 4. The relative error in the total surfactant mass versus time 
for mesh size /i = 2/64 = 0.03125. 
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Figure 19: Results for Example 4. The error \\{uB^h — i(o.5) (circles) and \\{us^h — 

'^S',2/i)||r/,(o.5) (stars) measured in the norm versus mesh size h. The time step size k = 0.625/i. 
The dashed line is proportional to h?. 



Figure 20: Results for Example 4. Condition number versus time for mesh size h — 2/64 = 
0.03125 and time step size k = h/8. 
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6 Conclusion 


We have presented a new finite element method for the coupled bulk-surface problem modeling 
the evolution of surfactants. The model is given by a time-dependent convection-diffusion equa¬ 
tion on the surface (or interface) coupled with a time-dependent convection-diffusion equation 
in the bulk. The concentration of surfactants in the bulk and the concentration of surfactants 
on the surface are coupled through a nonlinear term. For this problem we have proposed a 
space-time CutFEM. The variational formulation is written over the space-time domain. The 
space-time domain is divided into space-time slabs. Interpolation functions that are continu¬ 
ous in space but discontinuous from one space-time slab to another have been used and the 
discrete equations have been solved one space-time slab at a time. The space-time formulation 
allows easily for using unstructured meshes in time slabs which is useful when adaptive schemes 
are used. For multiphase flow problems adaptive schemes are of great interest because often 
quantities of interest are on the interface and errors are typically large close to the interface. 

We have added consistent stabilization terms in the weak formulation which guarantee that 
1) the system matrices have bounded condition number independently of the interface position 
relative to the background grid 2) the method is stable in case the problem is convection dom¬ 
inated. A great advantage with the presented method is that it is simple to implement both 
in two and three space dimensions since it relies on spatial discretizations of the geometry at 
quadrature points in time and thus much of the computations are the same as in the case of 
stationary problems. The numerical results demonstrate that the proposed method has optimal 
convergence order and conserves the total mass of surfactants independently of how the interface 
cuts through the fixed background mesh. 

In this paper, we have used a level set method and represented the interface with linear 
segments. The presented method can be used with any interface representation technique and 
it is a subject of future work to use a better approximation of the interface with the proposed 
method. Also, the velocity field has been given analytically in this paper. However, in future 
work we will couple the method to a flow solver solving for example the Stokes or the Navier- 
Stokes equations. 
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7 Appendix 

In some of the numerical examples the non-dimensional form of the equations have been used. 
Here we formulate the transport equations in non-dimensional form and report the changes in 
the weak formulation. 
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7.1 Non-dimensional form of the transport equations for soluble surfactants 

In the following we use the same non-dimensional variables as in [2]. Let L, /3^, and 
be the characteristic values for length, velocity, and surface and bulk surfactant concentration. 
The non-dimensional form of the surfactant concentration equations is given by 


dfUB + /3 • Vub — V • ~ ® 

in / X Di(t) 

(7.1) 

^ Da/coupling 

on r{t) 

(7.2) 

-'TidQ • ^'^ub = 0 

Pe 

on dD 

(7.3) 

dfUs + f3 ■ Vus + {divrf3)us - divp = /coupling 

on / X T{t) 

(7.4) 

with 

/coupling ~ O^Ub{^ Vs) Bi'U^' 


(7.5) 


Here the non-dimensional numbers Pe and Pcg are the bulk and surface (interfacial) Peclet 
numbers, defined with respect to the bulk and surface diffusivity ks and ks^ Da is the Damkohler 
number and Bi is the Biot number [2] 


Pe = 




Pes = 




Da = 


Uc 


Bi = 


kdL 


a = 




ks ^ ks ^ 

Due to the non-dimensionalization, the conserved total amount of surfactants is 


(7.6) 



usdv + Da 



Usds 


(7.7) 


7.2 Weak form 


We can write /coupling given in equation (7.5) in the form 


/coupling bsUB bgUs bsS^B'^S 


(7.8) 


with hs — hss = Oi and hs = Bi. Assuming that 65, hs^ and Da are positive constants we can 
multiply the bulk PDE ( |7.1[ ) with a test function ^ iL^(Di(t)) and the surface PDE ( |7.4[ ) 
as before with a test function hsvs G This yields after integration by parts and 

incorporating the boundary conditions the following weak form 


bB 

Da 


{dtUB,VB)Qi{t)+bs{dtUs, Vs)r{t)+(^(t, U, v)-(hBSUBUS, bBVB-bsVs)v{t) = 0 


with 


a(t, 7/, v) 


bB 

Da 


aB(t, ub,vb) + bsas(t, us, vs) + u, v) 


Mv = {vb, vs) G W 
(7.9) 

(7.10) 
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ass 9^8 before, as and as as before but with ks and ks replaced by Pe and Pcg, respectively, see 
equation (2.14). Note that compared to equation (2.12) there is only a change in the coefficient 
in front of as and (terms coming from the bulk PDE). The same changes are 

applied to the weak formulation of our space-time CutFEM described in SectionIn addition, 
the terms containing A and /i in equation (4.7) change since condition (2.10) changes to fTfl ). 
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